function data_set = counter_fact_mto(x, rs_series, type,part_eqm_flag)
    global bet0 bet1 tau a b w e  N_a N_w N_e w_max w_min alpha xi1 lb ub pi_e...
        sigmae elasts log_e_mu log_e_sigma skill_shock...
        base_rs shares wspread pe thetas pi_theta rho sig...
        unemp zeta  myslope nu xi eta gamma mu_a gamma b_init bet0_init bet1_init base_rs_init Pr share curve scale w_gamma pop_gr vouch_sh part_eqm

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Model Parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    b = x(1) ;
    xi =  abs(x(2)); 
    alpha = x(3);
    sig = abs(x(4));
    rho = abs(x(5));
    log_e_sigma = abs(x(6));
    gamma = x(14);
    tau = x(7);
    bet1 = (x(8));
    skill_shock = abs(x(9));
    pop_gr = 1.0893;
   
    % Preference Shocks
    theta_h = max(x(10), 1);
    theta_l = max(min(x(11), 1), 0); 

    % Preference Shocks Probility
    pi_theta = max(min(abs(x(12)), 1), 0);
    thetas = [theta_h, 1, theta_l];
    pi_theta = [1 - pi_theta, pi_theta];
     
    % Idiosyncratic shock variance
    sigmae = 1 /  abs(x(13))  ;  
    
    vouch_sh = 1; % MTO voucher cap
    part_eqm = part_eqm_flag; % Partial Equilibrium flag

    % These parameters can have arbitrary normalizations but aide
    % calibration
    base_rs = x(15); 
    bet0 = (x(16));
    scale = 1; 
    b = b*scale^(1+alpha);
    bet0 = bet0*scale^((1+alpha)/xi);
    bet1 = bet1*scale^((1+alpha)/xi);

    % Housing Parameters - Elasticities
    phi2 = x(end-1); 
    phi3 = x(end);
    phi1 = 1.75*3-phi2-phi3;
          
    % Optional parameters
    pe = 0.0; % Public education
    unemp = 0.0; % Unemployment benefits
    xi1 = 1; % Curvature on spillover
    w_gamma = 1; 

    % Housing Moments
    Hs_moments_base = [0.1934594	0.21735	0.2278978	0.2508063	0.3006039	0.2494808	0.2294	0.2151092	1.2530   1.2769		];    

    h_a = Hs_moments_base(1); 
    h_b = Hs_moments_base(5); 
    h_c = 1 - h_a - h_b;
    shares = [h_a; h_b; h_c];
       
    % Skill premium shock in SS
    eta = 1;

    N_at = 1; % Grid size for temporary component of ability
    N_ar = 40;
    N_a = N_ar * N_at;
    log_ar_sigma = sqrt(sig^2 / (1 - rho^2)); % Implied std. deviation
    log_ar_mu = 0.0;
    log_at_mu = 0.0;
    log_at_sigma = 0.0;
    
    % Parameters governing random component of wages
    N_w = 125; % Grid size for Income
    log_e_mu = 0.0;
    log_w_mu = 0.0;
    log_w_sigma = 0.9;
    corr_wa = 0.0;
    N_e = 120; % Grid over mean zero wage innovations
    
    lb = 0.0; % Lower bound of education
    ub = Inf; % Upper bound of education
   
    % Define convergence
    tol = 1e-4;
    max_It = 150;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Ability Process
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Grid and transition matrix for ability
    m_a = 4; % Std Deviations Away from SS
    [log_ar, Pr_ar] = tauchen(N_ar, log_ar_mu, rho, sig, m_a);
    ar = exp(log_ar);
    % Transitory component of ability
    at = 1.0; %[1.1, 0.9];
    pi_at = 1.0;

    % Combine into ability matrix
    a = kron(at, ar);
    a = a(:);
    log_a = log(a);
    log_a_mu = log_ar_mu + log_at_mu;
    % Std Deviation of log ability process
    log_a_sigma = sqrt(log_ar_sigma^2 + log_at_sigma^2);

    [~,a_ind] = sort(a);
    a = a(a_ind);

    Pr = repmat(kron(pi_at, Pr_ar), N_at, 1);
    Pr = Pr(a_ind, a_ind);
    norm_mean = dot((Pr'^50)*ones(N_a,1)/N_a, a);
      
    a = a / norm_mean ;
    base_rs = base_rs*scale;
    tau = tau*scale;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Wage Grid
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    high_low = 4.5;
    w = zeros(N_w,1);
    w(end) = high_low * log_w_sigma;
    w(1) = log_w_sigma*-high_low;
    zstep = (w(end) - w(1)) / (N_w-1);
    w(2:end-1) = w(1) + zstep * (1:N_w-2);
    wspread = exp(w);
    w = w + log_w_mu;
    w = exp(w);
    w_min = min(w);
    w_max = max(w);
    log_w = log(w);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Idiosyncratic shock of wages
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    e_min = logninv(0.00005, log_e_mu, log_e_sigma);
    e_max = logninv(0.99995, log_e_mu, log_e_sigma); 
    e = exp(log(e_min) + (log(e_max) - log(e_min)) * linspace(0, 1, N_e+1));
    e = e';
    pi_e = logncdf(e, log_e_mu, log_e_sigma);
    pi_mass = diff(pi_e); 
    pi_e = pi_mass / sum(pi_mass);
    e =  (e(1:end-1) + e(2:end))/2;
    e = e / dot(e, pi_e);   
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial Disribution of wealth and ability
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    mu = [log_a_mu, log_w_mu];

    % Covariance Matrix
    Sigma = [log_a_sigma^2 corr_wa*log_w_sigma*log_a_sigma;
             corr_wa*log_w_sigma*log_a_sigma log_w_sigma^2];

    [m1,m2] = meshgrid(log_a, log_w);
    Prob = zeros(size(m1));

    for i = 1:size(m1, 1)
        for j = 1:size(m1, 2)
            Prob(i,j) = mvnpdf([m1(i,j), m2(i,j)],mu,Sigma);
        end
    end

    Prob = Prob / sum(Prob(:));   
    
    % Elasticity
    elasts = [phi1, phi2, phi3 ]';
    Hs = ones(size(elasts));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial guess for rental rate and spillover
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Ss = [2.49569237238145;2.19773991908705;1.89032785122845]';
    Ss = Ss ./ dot(Ss, shares);
    rs = [base_rs*2, base_rs*1.5, base_rs]';
   
    tic()
    [Probss, Sss, Rss] = steady_state(Prob, Ss.^xi1, Hs, rs, tol, max_It);
    toc()

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % MTO
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    tic()
    data_set = generate_counterfactual_mto(Probss, Sss, Rss, Hs, tol, max_It, rs_series, type);
    toc()
    
    disp(["Completed MTO at scale" type])
end